To explore this we used the following parameters space :
constant_env=getAlllSummaries("unifPop50k")
I was not keeping track of the outpute rate so the convertion in real generation has to be done by hand:
constant_env$outputrate=10
constant_env=updateScale(constant_env)
Check the trajectories of the simulation for the population size for different mutation rate \(\mu\)(from left to right: \(\mu = 0, \mu = 0.001,\mu=0.01\)), different selective pressure \(\sigma_s\) (from top to bottomright: \(\sigma_s = 0.1, \sigma_smu = 0.2,\sigma_s=0.4,\sigma_s=1000\))and different carrying capacity \(K\)
plotAllTrajVar(constant_env,m=0.05,E=0,obs="N")
Check the variance at the
plotAllVariableSummaries(constant_env,E=0,estimate=eq2833b) #we use only simulations without noise
To explore this we used the following parameters space :
v=getAlllSummaries("growingDelta",exclude="outputrate")
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
noisy_env = rbind(v,constant_env[,colnames(v)])
noisy_env = noisy_env[ noisy_env$mu >0,]
To compare it we merge the two tables of results, and remove simualtion with no mutations
Check the trajectory off the effective population size:
plotAllTrajVar(noisy_env,m=0.05,E=0,obs="var_x")
Check the trajectory off the mean variance size:
plotAllTrajVar(noisy_env,m=0.05,E=0,obs="var_x",ylim=c(0,.005))
Check and compare the final variance :
plotAllVariableSummaries(noisy_env,E=0,estimate=eq2833b)
Some of the data seems to be missing because when \(\delta>\sigma\) we have extintions, has shown in the graph:
sigmas=unique(noisy_env$sigma)
cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
names(cols)=sigmas
tmpc=noisy_env[ noisy_env$sigma==0.4 & noisy_env$mu == 0.01 & noisy_env$m == 0.2 ,]
plot(tmpc$delta,tmpc$extinction,log="xy",type="n",ylab="extinction time",xlab=expression(delta),main=bquote(mu == 0.01 ~ omega == 0))
for(s in sigmas){
tmpc=noisy_env[ noisy_env$sigma==s & noisy_env$mu == 0.01 & noisy_env$m == 0.2 ,]
Ks=sort(unique(tmpc$K))
for(k in 1:length(Ks)){
t=tmpc[tmpc$K==Ks[k],]
r=tapply(t$extinction,t$delta,mean)
lines(sort(unique(t$delta)),r,type="b",pch=k,col=cols[as.character(s)])
}
}
legend("right",
legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
col=c(rep(1,length(Ks)),cols),
lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
pch=c(seq_along(Ks),rep(NA,length(sigmas)))
)
#dev.off()
Thus we introduce autocorrelation in the environment. To introduce autocorrelation we increase slightly \(omega\)
par(mfrow=c(1,3))
plot(1:500,environment(500,omega=1,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))
plot(1:500,environment(500,omega=2,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))
plot(1:500,environment(500,omega=3,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))
For this we re-run previous experiment but with \(\omega \in {1,2}\)
omega1=getAlllSummaries(c("omega1growingDelta"))
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
omega1=updateScale(omega1)
omega1$omega=1
omega2=getAlllSummaries(c("omega2growingDelta"))
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
omega2=updateScale(omega2)
omega2$omega=2
noisy_env_omega = rbind(v,constant_env[,colnames(v)],omega1[,colnames(v)],omega2[,colnames(v)])
Let’s look at the extinction again for the different omegas
par(mfrow=c(1,3))
omegas=unique(noisy_env_omega$omega)
sigmas=sort(unique(noisy_env_omega$sigma))
deltas=sort(unique(noisy_env_omega$delta))
Ks=sort(unique(noisy_env_omega$K))
for(o in omegas){
cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
names(cols)=sigmas
tmpc=noisy_env_omega[ noisy_env_omega$sigma==0.4 & noisy_env_omega$mu == 0.01 & noisy_env_omega$omega == o & noisy_env_omega$m == 0.2 ,]
xrange=range(noisy_env_omega$delta[noisy_env_omega$omega>0])
plot(tmpc$delta,tmpc$extinction,log="xy",type="n",ylab="extinction time",xlab=expression(delta),main=bquote(mu == 0.01 ~ omega == .(o)),xlim=xrange)
for(s in sigmas){
tmpc=noisy_env_omega[ noisy_env_omega$sigma==s & noisy_env_omega$mu == 0.01 & noisy_env_omega$m == 0.2 & noisy_env_omega$omega == o,]
for(k in 1:length(Ks)){
t=tmpc[tmpc$K==Ks[k],]
r=tapply(t$extinction,t$delta,mean)
lines(sort(unique(t$delta)),r,type="b",pch=k,col=cols[as.character(s)])
}
}
legend("bottomleft",
legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
col=c(rep(1,length(Ks)),cols),
lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
pch=c(seq_along(Ks),rep(NA,length(sigmas)))
)
}
We can thus reproduce he trajectories and variance equilibrium using simulation with \(\omega=1\) et \(\omega=2\)
Check the trajectory off the effective population size:
#pdf("trajectory_N_omega1.pdf",width=12,height=16)
plotAllTrajVar(omega1,m=0.05,E=0,obs="N")
#dev.off()
#pdf("trajectory_N_omega2.pdf",width=12,height=16)
plotAllTrajVar(omega2,m=0.05,E=0,obs="N")
#dev.off()
Check the trajectory off the mean variance size:
#pdf("trajectory_var_x_omega1.pdf",width=12,height=16)
plotAllTrajVar(omega1,m=0.05,E=0,obs="var_x",ylim=c(0,.002))
#dev.off()
#pdf("trajectory_var_x_omega2.pdf",width=12,height=16)
plotAllTrajVar(omega2,m=0.05,E=0,obs="var_x",ylim=c(0,.002))
#dev.off()
And check the variance at the end of the equilibrium:
constant_env$omega=0
omega1_wdz=rbind(omega1,constant_env[constant_env$mu>0,])
omega2_wdz=rbind(omega2,constant_env[constant_env$mu>0,])
omega1_wdz=omega1_wdz[omega1_wdz$sigma<10,]
omega2_wdz=omega2_wdz[omega2_wdz$sigma<10,]
plotAllVariableSummaries(omega1_wdz,E=0,estimate=eq2833b) #we use only simulations without noise
plotAllVariableSummaries(omega2_wdz,E=0,estimate=eq2833b) #we use only simulations without noise
par(mfrow=c(1,3))
vts=c(0.001,0.002,0.003)
cols=colorRampPalette(c("darkgreen","green"))(length(vts))
names(cols)=vts
omegas=1:3
for(o in omegas){
plot(1,1,xlim=c(1,500),ylim=c(-.1,1.5),type="n",xlab="t",ylab=expression(theta[t]),main=bquote(omega == .(o)))
for ( vt in vts)
lines(1:500,environment(500,omega=o,delta=.1,vt=vt),col=cols[as.character(vt)])
legend("topleft",legend=paste0("v=",vts),col=cols,lty=1)
}
moving_theta=getAlllSummaries("growingThetat")
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
moving_theta=updateScale(moving_theta)
par(mfrow=c(1,4))
rates=unique(moving_theta$vt)
sigmas=sort(unique(moving_theta$sigma))
deltas=sort(unique(moving_theta$delta))
Ks=sort(unique(moving_theta$K))
cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
names(cols)=sigmas
for(d in deltas){
tmpc=moving_theta[ moving_theta$sigma==0.4 & moving_theta$mu == 0.01 & moving_theta$delta == d & moving_theta$m == 0.2 ,]
#xrange=range(moving_theta$delta)
plot(tmpc$vt,tmpc$extinction,log="xy",type="n",ylab="extinction time",xlab=expression(delta),main=bquote(mu == 0.01 ~ omega == .(o)))
for(s in sigmas){
tmpc=moving_theta[ moving_theta$sigma==s & moving_theta$mu == 0.001 & moving_theta$m == 0.2 & moving_theta$deltas == d,]
for(k in 1:length(Ks)){
t=tmpc[tmpc$K==Ks[k],]
r=tapply(t$extinction,t$vt,mean)
lines(sort(unique(t$vt)),r,type="b",pch=k,col=cols[as.character(s)])
}
}
legend("bottomleft",
legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
col=c(rep(1,length(Ks)),cols),
lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
pch=c(seq_along(Ks),rep(NA,length(sigmas)))
)
}